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ABSTRACT 

We follow the collapse in axisymmetry of a uniformly rotating, supermassive 
star (SMS) to a supermassive black hole (SMBH) in full general relativity. The 
initial SMS of arbitrary mass M is marginally unstable to radial collapse and 
rotates at the mass-shedding limit. The collapse proceeds homologously early on 
and results in the appearance of an apparent horizon at the center. Although our 
integration terminates before final equilibrium is achieved, we determine that the 
final black hole will contain about 90% of the total mass of the system and have 
a spin parameter J/M 2 ~ 0.75. The remaining gas forms a rotating disk about 
the nascent hole. 

Subject headings: black hole physics - relativity - hydrodynamics - stars: rota- 
tion 



1. Introduction 

Recent observations provide increasingly strong evidence that supermassive black holes 
(SMBHs) of mass ~ 10 6 - 10 10 M o exist and that they are the central engines that power 
active galactic nuclei (AGNs) and quasars (Rees 1998, 2001). The dynamical formation of 
SMBHs, as well as the inspiral, collision and merger of binary SMBHs, are promising sources 
of long-wavelength gravitational waves for the proposed Laser Interferometer Space Antenna 
(LISA) (Thorne 1995; Schutz 2001). However, the actual scenario(s) by which SMBHs form 
is(are) still uncertain. Viable stellar dynamical and hydrodynamical routes leading to the 
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formation of SMBHs have been proposed (e.g., Begelman & Rees 1978, Rees 1984, Shapiro 
& Teukolsky 1985, Quinlan & Shapiro 1990, Rees 1998, 2001; Balberg & Shapiro 2002). In 
typical hydrodynamical scenarios, a supermassive gas cloud is build up from the multiple 
collisions of stars or small gas clouds in stellar clusters to form a supermassive star (SMS) 
(Begelman & Rees 1978). A supermassive gas cloud might be formed from scattered gas in 
the bulge that falls into the central region due to radiation drag (Umemura 2001). SMSs 
ultimately collapse to black holes following quasi-stationary cooling and contraction to the 
onset of radial instability (Zel'dovich & Novikov 1971; Shapiro & Teukolsky 1983). This 
scenario for forming a SMBH following SMS collapse was investigated in the 1960s and 
1970s, but the studies treated nonrotating configurations and assumed spherical symmetry, 
or employed approximate analytical models (see, e.g., Baumgarte & Shapiro 1999 for a review 
and references). 

However, SMSs are likely to be rapidly rotating (see, e.g., Loeb & Rasio 1994). Baum- 
garte & Shapiro (1999) recently performed a detailed numerical analysis of the structure 
and stability of a rapidly rotating SMS in equilibrium. Assuming the viscous or magnetic 
braking timescale for angular momentum transfer is shorter than the evolution timescale of 
a typical SMS (Zel'dovich & Novikov 1971; see New and Shapiro 2001 for an alternative), 
the star will settle into rigid rotation and evolve to the mass-shedding limit following cooling 
and contraction. (At mass-shedding, the angular velocity of gas at the equator equals the 
Kepler velocity). They found that all stars at the onset of quasi-radial collapse have an equa- 
torial radius R QAOGM/c 2 and a nondimensional spin parameter q = cJ/GM 2 « 0.97. 
Here J, M, c, and G are spin, gravitational mass, light velocity and gravitational constant. 
(Hereafter we adopt gravitational units and set c = G — 1). Because of the large value of q, 
it is uncertain whether the rotating SMS collapses directly to a black hole or forms a disk. 
Saijo et al. (2002) investigated the collapse of a rotating SMS in a post-Newtonian (PN) 
approximation, and concluded that effects of rotation do not halt the collapse and that a 
SMBH is likely to be formed. However, a PN calculation cannot follow collapse into the 
strong-field regime and cannot rigorously address the possibility of black hole formation and 
growth, and the final outcome. 

To clarify whether a SMBH forms as the endpoint of rapidly rotating SMS collapse and, 
if it does, to determine the final hole parameters, we performed a fully general relativistic 
numerical simulation of the collapse in axisymmetry. We demonstrate that the collapse 
indeed leads directly to a SMBH of moderately rapid rotation, with q ~ 0.75. 
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2. Computational Set-Up 

Simulations were performed using an axisymmetric code in full general relativity (Shi- 
bata 2000). This code was constructed from a 3D code (Shibata 1999b) using the so-called 
"cartoon method" in numerical relativity (Alcubierre et al. 2001): We solve the Einstein 
field equations in Cartesian coordinates with a grid of size (N,3,N) in (x,y,z), covering a 
computational domain < x, z < L and —Ay < y < Ay. Here N is a constant (3> 1), 
L the location of the outer boundary and Ay is the grid spacing in the ^/-direction, with 
an an axisymmetric boundary condition imposed at y — ±Ay; the spin axis is along z. 
We solve the hydrodynamic equations with an axisymmetric code in cylindrical coordinates. 
This hydrodynamical code was calibrated by numerous test calculations, including spherical 
collapse of dust, the stability of spherical stars, mode analysis of spherical stars, and the 
evolution of rotating stars (Shibata 1999b). We adopt maximal time slicing and approxi- 
mate minimal distortion as coordinate gauge conditions throughout the simulation (Shibata 
1999a, b). Formation of a black hole is determined by finding an apparent horizon. 

Violations of the Hamiltonian constraint and conservation of mass and angular mo- 
mentum are monitored as numerical accuracy check during the simulation. Total angular 
momentum J and total baryon rest-mass M* should be strictly conserved in axisymmetry. 
The gravitational mass M is not conserved, due to the emission of gravitational radiation, 
but the decrease in M is very small. 

Since the SMS spins up during the collapse, a nonaxisymmetric instability could arise. 
However, the PN study by Saijo et al. (2002) indicates that such an instability is not excited, 
at least when the equatorial radius of the collapsing star exceeds 10M. By restricting our 
analysis to axisymmetry, we can improve our resolution of the strong-field, central region 
where the black hole forms by a factor of > 10. 

To model the initial equilibrium SMS, we adopt a polytropic equation of state with the 
adiabatic index T = 4/3, setting P = Kp*/ 3 , where P and p are the pressure and rest-mass 
density. This prescription is appropriate for a radiation-dominated SMS equation of state. 
Here K is a constant whose value determines the mass; we scale out the mass by setting 
K — 1 (Baumgarte & Shapiro 1999). 

Although a spherical SMS with the adopted equation of state is unconditionally unstable 
to collapse, rotation can stabilize the star. We focus only on a rigidly rotating SMS at 
the mass-shedding limit. According to Baumgarte & Shapiro (1999), rotating SMSs with 
equatorial radii R < 640M are unstable against collapse. At R m 640M, the ratio of 
the rotational kinetic energy to the gravitational binding energy T/W is about 0.009 and 
q « 0.97. 
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We chose a SMS with R « 620M, a star which is located just beyond the critical point 
of instability. In this case, we have T/W ~ 0.0088 and q ~ 0.96, values which are nearly 
equal to those of the critical configuration (see Table 1). A second simulation was performed 
for a more compact star with R R3 570M, but our results were essentially unchanged. To 
accelerate the collapse, we depleted the pressure by 1% initially. During the evolution, we 
adopted a T-law equations of state according to P = (T — l)pe where e denotes the specific 
internal energy of the fluid, thereby treating the gas as an ideal, adiabatic fluid in which 
cooling can be ignored on a dynamical timescale (Linke et al. 2001). 

Since the equatorial radius decreases by a factor of ~ 1000 (from ~ 600M to < M), 
using a fixed uniform grid with sufficient resolution for all epochs would be computationally 
prohibitive. Instead, we adopted a uniform grid with decreasing grid spacing and increasing 
grid number as the collapse proceeded in order to guarantee adequate resolution up to the 
formation of an apparent horizon. For the early stages, where the radius at the equatorial 
surface R m exceeds 150M, we used iV = 300 grid points in the x and z directions. Here, we 
identify R m as the radius at which p = 10 _6 p max , where p max is the maximum interior density. 
Initially, the grid covers the equator with about 200 grid points. The outer boundaries along 
the x and z axes are located at L m 930M. 

Since T/W is small and T = 4/3, the collapse proceeds in a homologous manner in 
the central regions during the early stages (Shapiro & Teukolsky 1979; Goldreich & Weber 
1980; Saijo et al. 2002). Taking advantage of this behavior, we rezoned by moving the outer 
boundary inward and decreasing the grid spacing, keeping N(= 300) fixed. All quantities 
in the new grid are calculated using cubic interpolation. We discarded the outermost com- 
putational domain, but the discarded baryon rest-mass is very small (less than 10~ 4 of the 
total) when R m ^ 150M. We repeated this procedure twice until homology breaks down at 
R m <> 150M. After this stage, the collapse timescale in the central region is much shorter 
than in the envelope. Consequently, we increased the grid number N and decreased A, 
monitoring the lapse function at the center otQ. Specifically, we set N and L as follows: for 
a > 0.90, we set N = 600 and L « 233M; for 0.70 < a < 0.90, we set N = 900 and 
L « 155M; and for 0.30 < a < 0.70, we set N = 1200 and L « 116M. With this treatment, 
the discarded fraction of the baryon rest-mass is only ~ 0.15% down to a = 0.3. 

For ao < 0.3, the central density profile becomes very steep. To see its dependence 
on grid resolution during black hole formation, we carried out simulations using different 
combinations of iV and L to refine the grid, with the restriction that N < 2400. The set-up 
of each simulation was as follows: (A) for a < 0.30, we took iV = 1200 and L 116M, with 
no regridding at a = 0.3; (B) for a < 0.30, we took iV = 1200 and L 58M, regridding 
once at a = 0.3; (C) for 0.05 < a < 0.30, we took N = 1200 and L w 58M and for 
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a < 0.05, we took N = 1200 and L w 29M; (D) for 0.05 < a < 0.30, we took N = 1200 
and L pa 58M and for a < 0.05, we took iV = 2400 and L pa 58M. For cases (C) and 
(D), we carried out regridding at «o — 0.3 and 0.05. The minimum grid spacings for cases 
(A)-(D) were 0.097M, 0.048M, 0.024M and 0.024M, respectively. The results for case (D) 
are the most reliable, but the results for the four cases do not differ significantly. 

For cases (B)-(D), the outer boundaries reside deep inside the stellar surface. Hence, 
the fluid at large equatorial radius is discarded. When we allow L pa 58M (29M), the loss 
of total baryon rest-mass becomes ~ 2% (4%). Since the outer envelope has larger specific 
angular momentum, setting L pa 58M (29M) implies that ~ 8% (18%) of the total angular 
momentum is lost. As shown below, however, discarding some mass in the outer region is a 
tiny effect on the formation and evolution of a SMBH in the central region. 

3. Numerical Results 

In Fig. 1, we display snapshots of density contours and velocity vectors in the x - z 
plane at selected times. Here, the velocity field is defined as where is the four 

velocity. The collapse proceeds nearly homologously in the central region during the early 
stages. However, for R m <> 150M, the effects of rotation and general relativity modify this 
property, when the collapse near the central region accelerates significantly. 

In Fig. 2, we show the time evolution of the central conformal factor (tpo = [det^j)] 1 / 12 
at r = where 7^ denotes the three-dimensional spatial metric) and «o for cases (A)-(D). 
We find that the collapse proceeds in a runaway manner in the final stages, although the 
time development depends somewhat on our adopted choice of time slicing. 

Since ip 2 ■ A measures the proper physical length of the grid spacing, maintaining the 
resolution requires changing A as A oc ip~ 2 (i.e, N oc ip 2 ). However, ip diverges very sharply 
at late times, so increasing A by a factor of a few does not improve the resolution much. 
Accordingly, we terminated the simulation when the conservation of M was violated by 
~ 10% at t/M pa 30636. 

There are two reasons for this runaway behavior at the center. One is that the equation 
of state is very soft, which produces stars with very centrally condensed structures. The 
collapse timescale in the central region is then much shorter than that in the outer region. 
The other reason is our choice of a coordinate (shift) condition, for which the resolution 
near the center deteriorates ( "grid sucking"; see Shibata 1999b). Integrating to a final 
equilibrium state evidently requires different time and/or spatial gauge conditions. 
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We find an apparent horizon forms for t/M ^ 30630. In Fig. 3, we show the mass of 



A denotes the area of the apparent horizon (e.g., Cook & York, 1990). In addition to Mah, 
we plot the total baryon rest-mass inside the apparent horizon M*ah for case (D) (dotted 
curve), which agrees approximately with Mah- Figure 3 indicates that at the end of the 
simulation about 60 % of the total rest-mass already has been swallowed into the black hole. 
Clearly, collapse of a rapidly and uniformly rotating SMS leads to formation of a SMBH. 
However, our final snapshot is not the final state, because Mah is still increasing. 

It is possible to estimate what the final mass and spin of the black hole will be once 
it settles into equilibrium. Define the specific angular momentum according to j = hu v , 
where h(= 1 + e + P/p) is the specific enthalpy. In an axisymmetric system, the integrated 
baryon rest-mass of all fluid elements with j less than a given value jo, m *(jo) (the angular 
momentum spectrum), is conserved in the absence of viscosity (Stark & Piran 1987): 



Consider the innermost stable circular orbit (ISCO) around the growing black hole at the 
center. If j of a fluid element is less than the value at the ISCO (jisco), the element will fall 
into the black hole eventually. Now the possibility exists that some fluid can be captured 
even for j > jisco, if it is in a noncircular orbit. Ignoring these trajectories yields the 
minimum amount of mass that will fall into the black hole at each moment. The value 
of jisco changes as the black hole grows. If jisco increases, additional mass will fall into 
the black hole. However, if jisco decreases, ambient fluid will no longer be captured. This 
expectation implies that when jisco reaches a maximum value, the dynamical growth of the 
black hole will terminate. 

To analyze the growth of the black hole mass, we generate Figs. 4 and 5. In Fig. 4, 
we show the angular momentum spectrum at t — 0. To verify that the spectrum is well 
preserved during the simulation, we also plot the spectrum at the time when the apparent 
horizon is first formed. In Fig. 5 (a), we plot q* = J {j) / ™*{j) 2 as a function of m*(j)/M*. 
Here, J(j) denotes the total angular momentum with the specific angular momentum j less 
than a given value jo and is defined according to 



the apparent horizon as a function of time. The mass is defined as Mah = 



a/ A/lQn where 




(1) 




(2) 



Now, J(jisco)/' m *(jisco) 2 and m^(jisco) m ay be approximately regarded as the instanta- 
neous spin parameter and mass of a black hole. Therefore, the dotted curve in Fig. 5(a), 
derivable from the initial stellar profile, may be interpreted as an approximate evolutionary 
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track for the angular momentum parameter of the growing black hole. The numerical results 
for various resolutions indicate that this is close to the actual track followed by the hole. 
Thus, we can assume that our assumptions made in this analysis are adequate. The solid 
curve in Fig. 5(a) shows that with the increase of the black hole mass, the spin parameter 
also increases. 

If we assume that m*(j) and g* are the mass and spin parameter of the black hole 
and that the spacetime can be approximated instantaneously by a Kerr metric, we can 
compute jisco (see, e -g-> chapter 12 of Shapiro & Teukolsky, 1983). In Fig. 5(b), we show 
Jisco[^*(j), Q*(J)] as a function of m*(j). For m*(j) < 0.9M*, jisco is an increasing function 
of the mass. This implies that the mass of the black hole should increase to ~ 0.9M*. 
However, at m^j^/M* 0.9, jisco reaches a maximum, as jisco/ m *(j) steeply decreases 
above this mass fraction. Thus, once the black hole reaches this point, it will stop growing 
dynamically. Figure 5(a) shows that at this stage, ~ 0.75. Therefore, at the end of this 
collapse, (1) about 90% of the total mass will form a SMBH and (2) the spin parameter of 
the Kerr hole at the end of the collapse will be ~ 0.75. 



4. Summary 

We performed a fully relativistic numerical simulation in axisymmetry of the collapse 
of a uniformly rotating, marginally unstable SMS. Our simulation terminates when roughly 
60% of the mass has been swallowed by the SMBH. We estimate that about 90% of the total 
mass of the system will be consumed by the end of the collapse. The spin parameter q of 
the final Kerr SMBH is likely to be ~ 0.75 at the end of the dynamical collapse phase. Most 
of the remaining gas will reside in an ambient disk about the central hole. 

To follow black hole growth to a final equilibrium state will require different coordinate 
gauge choices for the lapse and shift functions. It may also require the use of an "horizon 
excision" boundary condition (Unruh, unpublished; Thornburg 1987; Seidel & Suen 1992). 
To simultaneously follow the extended disk may necessitate employing a nested grid or 
adaptive mesh refinement. These issues are ripe for future exploration. 

We are grateful to T. Baumgarte, Y. Eriguchi, H. Shinkai, H. Susa, M. Umemura, and 
K. Uryu for discussion. Numerical computation was performed on FACOM VPP5000 in the 
data processing center of National Astronomical Observatory of Japan. This work was in 
part supported by a Japanese Monbu-Kagaku-sho Grant (No. 13740143) and by NSF Grant 
PHY-0090310 and NASA Grants NAG5-8418 and NAG5-10781 at the University of Illinois 
at Urbana-Champaign. 
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Table 1: Parameters of the initial SMS. 



R/M 


T/W 


q M 


M* 


n 


Pc 


initial data 622 


0.0088 


0.96 4.566 


4.566 


1.40e-5 


7.84e-9 



Note. — From left, the radius at the equator, ratio of the kinetic energy to the gravitational binding 
energy, gravitational mass, baryon rest-mass, angular velocity and central density. All the quantities are 
shown in units of c = G = K = 1. 
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Fig. 1. — Snapshots of density contours and velocity vectors in the x-z plane at selected 
times for case (D). The contours are drawn for p/p max = lCT ' 4 -' (j — ~ 15), where p max 
denotes the maximum density at each time. The fourth figure is a blow-up of the third one 
in the central region: The thick solid curve at r 0.3M shows the location of the apparent 
horizon. 
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Fig. 2. — The central value of the conformal factor ipo and lapse function a>o as a function of 
time. The solid, dotted, dashed, and dotted-dashed curves denote the results for cases (D), 
(A), (B) and (C). The results for (C) and (D) are not distinguishable. 
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Fig. 3. — Evolution of the mass of apparent horizon (solid curve) and baryon rest-mass 
inside the apparent horizon (M*ah, dotted curve) as a function of time for case (D). Mass 
and time are shown in units of the total gravitational mass M. For comparison, we also plot 
the mass of apparent horizon as a function of time for cases (A) (crosses), (B) (triangles) 
and (C) (circles). The results for cases (B), (C) and (D) essentially agree. 
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Fig. 4. — The specific angular momentum spectrum at t « (dotted curve) and at the first 
formation of an apparent horizon at t « 30630M for cases (D) (squares), (A) (crosses) and 
(B) (triangles). 
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Fig. 5. — (a): J(j)/m*(j) 2 as a function of m*(j) (dotted curve). The squares, circles and 
triangles show Jah/M^ ah and M*ah at select times for cases (D), (C) and (B), respectively. 
Here, Jah and M*ah are the angular momentum and baryon rest-mass inside the apparent 
horizon, (b): jisco/M as a function of m*(j)/M*. The cross marks the maximum. 



